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Abstract 

We present and analyze an approximation scheme for the two-dimensional 
game p-Laplacian in the framework of viscosity solutions. The approximation 
is based on a semi-Lagrangian scheme which exploits the idea of p-averages. 
We study the properties of the scheme and prove that it converges, in par- 
ticular cases, to the viscosity solution of the game p-Laplacian. We also 
present a numerical implementation of the scheme for different values of p; 
the numerical tests show that the scheme is accurate. 

Keywords: p-Laplacian, Tug of war game, Hamilton- Jacobi equations, 
semi-Lagrangian scheme, convergence, viscosity solutions. 



1. Introduction 

The game p-Laplace operator has been recently introduced in [25] to 
model a stochastic game called tug-of-war with noise. Part of the interest 
for this class of operators arises from the fact that it includes, as particular 
cases, the operator in the Aronsson equation [I], the infinity Laplacian |24j, 
the motion by mean curvature operator [XT] , and for p = 2, a multiple of 
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the ordinary Laplacian. For the connections between these operators and 
differential games see also [ID] . 

The equation associated to the game p-Laplacian has the same solutions 
as the variational p-Laplacian only in the homogeneous case, and has the 
advantage in the non-homogeneous case of being a combination of other p- 
Laplacians. In particular, we would like to stress the fact that even for the 
non-homogeneous case the game oo-Laplacian is the limit as p — > oo of the 
game p-Laplacian. 

Our work strongly relies on the general philosophy illustrated in the paper 
of Peres and Sheffield [25] , which indicates that to study p-harmonic functions 
one can look at discrete versions of the p-operators; these correspond to 
stochastic processes with paths that are nonlinearly averaged, ranging from 
motion by mean curvature to Brownian motion to diffusions generated by 
the Arronsson operator. It is important to note that our work is in the 
framework of weak solutions in the viscosity sense (see [9J for an introduction 
and [7] for a guide to viscosity solutions for second order problems) . However, 
the difference between the analytic case of Peres and Sheffield [25], and our 
approximation scheme is in the fact that we need the value at a fixed point 
to depend only on a discrete number of values in space. In this respect, our 
construction starts from an approximation of p-averages in order to keep a 
strong link between the continuous and the discrete operator. Specifically, 
one of our goals is to prove that our scheme is consistent with the continuous 
operator proposed by Peres and Sheffield. 

Semi-Lagrangian schemes for nonlinear Hamilton- Jacobi equations have 
been studied and analyzed by several authors. The starting point is the dis- 
crete version of the characteristic method which leads to a discrete Lax-Hopf 
formula for first order Hamilton- Jacobi equation [IT] . It is interesting to note 
that semi-Lagrangian schemes allow for large time steps still satisfying stabil- 
ity conditions. A comprehensive introduction to semi-Lagrangian methods 
for linear and first order Hamilton- Jacobi equations is contained in the book 
by Falcone and Ferretti [12]. For second order problems some results have 
been obtained for stationary and evolutive equations related to stochastic 
control problems in [5], whereas a presentation of the treatment of second 
order terms in SL schemes has been given in [TJ] . Recent extensions to mean 
curvature driven flows have been analyzed in [6j. 

In the homogeneous case the game p-Laplacian coincides with the varia- 
tional p-Laplacian for which several approximation methods have been pro- 
posed. Some of these schemes are based on finite elements and show conver- 
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gence, also establishing a priori error estimates, see e.g. the paper by Barrett 
and Liu [3|. However, finite elements are not the most popular techniques 
for nonlinear degenerate equation. Finite difference approximation schemes 
for degenerate second order equations have been proposed and analyzed by 
Crandall and Lions in |8j, and in several papers by Oberman [2"T| [221 123] ■ 
More recently, finite volumes schemes have been presented by Andreyanov, 
Boyer and Hubert in p], for the variational p-Laplacian. 

In the non-homogeneous case the game p-Laplacian interprets the non- 
homogeneity, /, as a multiple of a running payoff for one of the players in 
a two-player, zero-sum game while for the variational p-Laplacian the non- 
homogeneity is interpreted as a potential. From the numerical point of view, 
the variational p-Laplacian is typically studied with homogeneous Dirichlet 
boundary conditions. In this case, one has Poincare inequalities when the 
region and solutions are smooth, which help in the proofs of convergence and 
error estimates. The game p-Laplacian has fewer tools, relying essentially 
on monotonicity properties and on the notion of viscosity solutions. For this 
reason we believe that our results will be a useful contribution to the theory. 

The paper is organized as follows. In Section 2 we introduce the formal 
definition of the game p-Laplacian, and give a precise definition of viscosity 
solutions for our problem, as well as a brief description of the stochastic game 
tug-of-war with noise. We formulate our approximation scheme in Section 3, 
where we also give insights into its construction. The analysis of the prop- 
erties of the scheme, and a proof of convergence for the homogeneous case 
when p > 2 are presented in Section 4. In Section 5 we discuss in detail the 
numerical implementation of the scheme on a rectangular grid, and we con- 
clude with some numerical tests in Section 6. For the sake of completeness, 
we also add in Section 7 a technical appendix on some elementary properties 
of the p-average of finite sets of real numbers. 

2. Game p-Laplacian 

The p-Laplace operator, which we refer to as the variational p-Laplacian, 
for 1 < p < oo, is defined by 




(1) 



whereas for p = oo, traditionally is given by 

du du d 2 u 
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The subject of our numerical study is the Dirichlet boundary value prob- 
lem for the so-called game p-Laplacian introduced by Peres and Sheffield [25J , 
which for 1 < p < oo reads as follows: 



A°u = f in 



p 



u = F on d£l, 



(2) 



where 

ApU := - |Vn| 2 - p div (\Vu\ p - 2 Vu). (3) 

We require / and F to be continuous in their domain of definition. Ad- 
ditionally, we assume / either identically equal to zero or never zero. In the 
sequel, we will then always consider the cases / = or / > ( without loss 
of generality). Here Q C 1R™ is a bounded smooth domain. 

Note that for p = 2, one has AJf = 5A2, that is one-half of the Laplacian, 
which is the infinitesimal generator for a Brownian motion. 

If u is a smooth function, by expanding the derivatives, we obtain 

A° ^ A + ^ ~ ^ IV |~ 2 ^ U ^ U (4) 
p p p ^ri dxi dxj dxidxj ' 

therefore, by taking the limit for p — > oo, one is naturally lead to the following 
definition of the game oo— Laplacian: 

^ G ,^ |_ 2 du du d 2 u 
00 ^ dxi dxj dxidxj 

The game 1-Laplacian is defined in terms of the Laplacian and the game 
oo— Laplacian: 

A?u := A 2 u - Agu. (6) 

We would like to point out that with this notation, the expansion in Q 
allows us to think of A^ as the convex combination of the two limiting cases, 
that is 

= - Af + - A^, (7) 
F p q 

with q the conjugate exponent of p (i.e. - + ^ = 1). 

At the points where Vu ^ 0, the game 1-Laplacian and the game oo- 
Laplacian can be thought as the second derivative in the orthogonal direction 
of Vu and in the direction of Vu, respectively. That is, 

Af u = \Vu x \~ 2 < D 2 uVu L , Vu 1 - >, (8) 
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and 

A^u = \Vu\~ 2 < D 2 u Vu, Vu >, (9) 

where D 2 u denotes the Hessian matrix. 

The variational p-Laplacian can be obtained as the Euler-Lagrange equa- 
tion of an energy functional, a fact that does not hold for the game p- 
Laplacian. Additionally, while the variational p-Laplacian is degenerate el- 
liptic for 2 < p < oo and singular for 1 < p < 2, the game p-Laplacian 
is singular for every p ^ 2, so suitable definitions of viscosity solutions are 
needed. 

Juutinen and al. [IB] have shown that for the variational p-Laplacian, 
when 1 < p < oo, the notions of viscosity solution and weak solution are 
equivalent. The interested reader can find in the survey [7] a number of 
results on viscosity solutions for second order problems. Note that, in the 
homogeneous case, i.e. when / = 0, the solutions of the two operators agree 
with each other. 

Various definitions of viscosity solutions for the game p-Laplacian can be 
given and are found in the literature. The most suitable for our treatment is 
the one obtained by following the definition in the classical paper of Barles 
and Souganidis [2]. 

In what follows, we will restrict ourselves to the two-dimensional case, 
that is we will take ft C M 2 . 

Definition 2.1. Consider a smooth domain ft C M. 2 , and let 1 < p < oo, q 
such that 1/p + l/q — 1. If f is a continuous function, we say that an upper 
semi- continuous function [respectively, lower semi- continuous] u : Q — )■ K is 
a viscosity subsolution [supers olution] of 

- A^u(x) = f(x) in ft, (10) 

if for any G C 2 (ft) such that u — <fi has a local maximum [local minimum] 
at x G ft, we have 

(i) -A^(x) < f(x) [-Ap(f)(x) > f(x)} if V0(x) ^ ; 

(ii) if Ai < A2 denote the eigenvalues of D 2 (p(x), then: 

_Ai _ A2 ^ Ai _ A2 > fix)] ifV<f>(x) = and p> 2; 

p q q p 

-—- — <f(x) [-—- — > fix)] ifVMx) = andl<p<2. 
q p p q 
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Remark 2.1. Part (ii) of the definition of viscosity subsolution [supersolu- 
tion] is implied by the condition: 



-Af<f>{x) < f{x) [-Af</>(x) > f{x)} whenever V0(x) 
This is a consequence of the fact that 

Ai A2 



0. 



and 



Ai At 

i £ < 

p (1 

Ai A2 



x < 
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P 



< -A 



G , 



Q 



V 



(x) < - 



Ai A2 
P Q 



ifp>2, 



if 1 < p < 2. 



Uniqueness for viscosity solutions of nonlinear operators that are singular 
at isolated points, typically does not depend on the particular value one 
assigns to these points as long as this is chosen in a consistent manner (see for 
example Section 9 in [7]), additionally our numerical results show numerical 
convergence to solutions that verify (ii)'. Therefore, we will use the following 
definition for viscosity solution of the game p-Laplacian: 

Definition 2.2. A function u is a viscosity solution of (T(fy , for 1 < p < 00 



if u is a subsolution and a supersolution according to (i) of Definition 2.1 



and (ii)' in Remark 2.1 



As we said in the introduction, the main ingredient of our approximation 
is the way we discretize using p-averages. Let us recall the notion of p-average 
of a set of numbers. 

Definition 2.3. Given a finite set of real numbers, S = {s\, S2, , s m }, we 

denote by A P (S) the p-average of its elements, that is A P (S) is such that 



i=i 



if 1 < p < 00, 



and 



A 1 (S) 



max Sj + min s,- 

Sj&S J Sj&S J 



median of S. 



(11) 
(12) 
(13) 
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Since the median of an even number of points is not uniquely defined, 
in (13) we follow tradition and take it to be the average of the two middle 
points. 

Note that by convexity A P (S) above is unique for 1 < p < oo. For p = 2, 
A2(S) is the arithmetic mean of the numbers in the set S: 



1 m 



m 

We end this section with a brief description of the two-player, zero-sum 
game called tug-of-war with noise. In this game, we fix a parameter e > 0, 
a 1 < p < oo, a domain Q, a continuous running cost h : Q — > R with 
h > 0, and a continuous exit cost F : dVt — > K, as well as an initial position 
x = x G Q. A token is placed at x , and at each stage, k, a fair coin is flipped. 
The winner of the coin flip picks any direction vector, v^, with \vk\ < e, to 
which a random noise vector Zk is added. The vector Zk is equally likely to 

be one of the two vectors orthogonal to Vk of length \ |t>J. Then the 

V p-i 

token is moved to Xk = Xk-i + Vk + Zk and the play continues until the token 



is within a distance e*= [ 1 + \j e from dVt. In that case, the winner 

P-IJ 

picks Xk G dVL with \xk — Xk-x\ < e". The payoff to one of the players, say 

fe-i 

player I, from the other player, say player II, is F{xk) + e 2 h(xi). Under 

i=0 

various conditions, Peres and Sheffield [25] show that when both players 
choose optimal strategies, as e — > there is an expected value u(x) = u(x ) 
which solves the boundary value problem for the game p-Laplacian equation 
1 2 

given in (2) with / = -h. Note that as p — > oo the noise vector disappears, 
U 1 

while for p — 2 it has the same length as the chosen direction vector, resulting 
in a two-dimensional random walk. 



3. Construction of the approximation scheme 

We arrive at our approximation scheme inspired by the work of Peres and 
Sheffield [25], where the game interpretation of the game p-Laplace operator 
is based on averaging over non-Markovian paths, by asking the question if 
the p— game operators have an averaging characteristic, and if this can be 
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captured by some quantity. Having this in mind, looking at the numerical 
approximations of Oberman [21, 22J of the 1-Laplacian and the oo-Laplacian, 
as well as at the standard central difference approximation of the 2-Laplacian, 
we notice that they all can be rewritten in terms of their corresponding p- 
average. 

From these observations, we arrive at the conclusion that there should 
be an inherent averaging characteristic in all the p-operators, and that the 
notion of p-average is a possible candidate for the correct quantity describing 
it. We refer to the work [15] for an analytical result on averaging properties 
of the p-Laplacian in terms of a continuous p-average (see also |2U]). 

To stress how the notion of p-average comes into play quite naturally, 
when dealing with approximation schemes for the game p-Laplacian, let's 
look at some known approximation schemes and rethink them in terms of 
p-averages. To avoid cumbersome notations we present the approximation 
scheme in M. 2 but the extension to the general multidimensional case follows 
along the same lines. 

For p = 2, using standard central differences, at a point x = (x\,x 2 ), for 
/iGl small, we have 

u(xx, x 2 ) « M x i + K x 2 ) + u(x x , x 2 + h) 
+u(x 1 - h, x 2 ) + u(xi, x 2 — h) — 4 u(xi, x 2 )) , 

which, reordering terms in a suitable way, gives 

A««(x)«-^[A 2 (C h (x,«))-«(x)]. 

In general, we will denote by C/j(x, u) the set of values used to compute the 
approximation at a point x, in this case C/j(x, u) = {u{x\ + h, x 2 ),u(xi,x 2 + 
h), u(xi — h, x 2 ), u{x\, x 2 — h)}. 

For p = oo, the scheme in [21] can be rewritten as 

~ \ 2 [Ax,(C h (x,u)) - u(x)] , 

where now C/i(x, u) is a discrete set of values of u on the sphere of radius 
h centered at x, and the distribution and number of points on the sphere 
influences the accuracy of the approximation in a fundamental way. 
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It is relevant to mention that also for p = 2 one could pick as C^(x, u) 
a larger set of values of u on the sphere, but since the Laplacian is a linear 
operator this would not increase the accuracy. 

A similar scheme that uses the 1— average can be constructed in the case 
p — 1, in view of the interpretation as second directional derivative given by 
(|8]), see [22J for the parabolic case. 

The generalization to the game p-Laplacian of these interpretations using 
averages suggests the following approximation: 



«(x) « ^ [MCh(x, «)) - u(x)] , (14) 



where again C^(x, -u) would be a suitable discrete set of values of it on the 
sphere of radius h centered at x. 

We are then lead to the following approximation scheme for the Dirichlet 
problem Q: 

S(p, x, it(x), u) — in f2, (15) 

where the positive discretization parameters are represented by the vector 
p := (h,A9) (with h the spatial step and AO the angular resolution), and 
S : [0, 1) x (0, 7r/2] xflxRx L°°(n) — ► E is defined as 

f [A p (C(x,«;a)) - u(x)] - /(x) in fi, 

S(p,x.,u(x),u) = <j (16) 
[ w(x) - F(x) on 

Here, if cfo < oo denotes the diameter of Q, a = a(x) is a dilation parameter 
such that < a(x) < dist(x.,dQ) < dn- Finally, C^ e (x,u;a) is now a 
suitably chosen discrete set of values of u, taken on the sphere of center x 
and radius ha, associated to the angular resolution A#. 

There is some freedom in how to choose the points in Cfr 9 (x.,u; a), but 
as in [13] we follow a standard discretization of the sphere of radius ha 
centered at x, and take yj = x + h a r j (with A6 = it /(2 m) and r j = 
(cos iA8, siniA0)), so that 

Q Ae (x, M ; a) = {u(y t ), i = 0, ...,4m - 1}. (17) 

Note that with this choice of m, if a direction r-j is in the set of admissible 
directions so is its opposite, — r i; as well as its orthogonal and its reflections 
with respect to each of the axes. Also, note that with our choice of a we 
have that G for every i, so our set C^ 9 {x, u; a) is well-defined. 
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4. Study of convergence 



Let us analyze the convergence of our approximation scheme using the 
framework provided by the classical result of convergence for fully nonlinear 
second order elliptic equations of Barles and Souganidis presented in [2] . As- 
suming that a comparison principle is available for the exact solution, in this 
approach convergence to viscosity solutions is implied by the monotonicity, 
stability and consistency of the scheme. 

We prove monotonicity for the general case, consistency for p > 2, and 
stability for the case / = 0. Therefore, we have formal convergence of the 
scheme for the case p > 2 and / = 0. Nevertheless, the numerical experiments 
we run, and which we illustrate in Section [6j show convergence in the general 
case. 

We follow and define G p : T 2 x R 2 x E x -> E, where T 2 is the 
set of 2 x 2 real symmetric matrices, by first introducing the function H p : 
T 2 xK 2 xfl4ft 



# p (M,r,x) = 
and then setting 



1 r r I r r 

-<M—,— > — <M— ,— >-/(x) if r ^0 
p | r- 1 - 1 | r- 1 - 1 q |r| |r| 

-itr(M) - /(x) ifr = 0, 



if p (M,r,x) for x G fi, 
G p (M,r, M ,x) = <( ' (18) 
w(x) - F(x) for x G cKl 

In this notation, the Dirichlet problem ^ is expressed as 

G p (D 2 u, Du, u, x) = for x e H. (19) 



Remark 4.1. A viscosity solution of (19) is a function u that verifies Defi- 



nition 2.2 in Q and satisfies the boundary conditions 



m&x(H p (D 2 u, Du, x), u — F) > on <9f2, 
mm(H p (D 2 u, Du, x), k — F) < on dVL. 

The ellipticity of G p is a trivial consequence of its definition: 
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Lemma 4.1. G p is elliptic, namely for all (r,n,x) 6 I 2 x R x f2 we have 

G p (M,r,«,x) < G p (iV,r,n,x), 
/or a// M, NeT 2 such that M — N is positive semidefinite. 

Monotonicity of the scheme stems from the monotonicity property of the 



p-average, which we derive, for the convenience of the reader, in Lemma 
of the Appendix. 

Theorem 4.2. Let u,v E L°°(tt), if n(x) > u(x) in then for all p > I, 
p G [0, 1) x (0, 7r/2], xgH, and t G R it holds 

S(p,x.,t,u) < S(p,x,t,v). 

Proof. If x e <9f2 then S(p, x, £, u) = t — -F(x) = S^p, x, £, n) and the claim is 
clearly true. If x e ^ we have that 

S(p,x,t,u) = [A p (C ft A9 (x,«;a)) - f] - /(x) 

^ [A>(CfW;")) - t] - /(x) = S(p,x,M), 



a- 



since A p (C^ 6> (x, m; a) ) > A p (C^ e (x, u; a)), thanks to the assumption u>v 
in VL and Lemma u.2l □ 

To prove consistency, we start by showing that in the case p > 2 our 
approximation has the correct behavior in the internal points of the domain. 

Theorem 4.3. Let p > 2. For all x e and e C°°(fi) ; we aawe £aa£ 

f A«0(x) z/V0(x)^O 

A^( X ) if V0(x) = 



Proof. Assume V0(x) ^ 0, and denote by ei = (1, 0); without loss of gener- 
ality, we can assume x = (0,0) and V0(x) = |V0(O)| ei. Equation ^ then 
gives 

A^>(0)=cW(0), (21) 

while (J8f yields 

Af0(O) = <9 22 0(O). (22) 
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The M:=4m points in the set Cf^ (x, 0; a) are now on a sphere of center 
and radius ah, with < a < dn, see (17), that is they are given by 

£ = a/i (cos 0,-, sin Bj) for j = 1..M, (23) 

where the 0j are uniformly distributed angles verifying \9j+i — 0j \ = A9. We 
use Taylor's expansion to obtain 

</>(£;) = 0(0) + V0(O) • & + i < £>V(0)£;,& > +o(a 2 /i 2 ), 

and by Lemmas |7.3 and |7.1| it is enough to show that the p-average of 
-^-^ (vcf>(0) ■ & + i < D 2 (f)(0)^, tends to A^0(O) as /i and A0 tend 

to 0. 

If p — oo, consistency is proven as in Oberman [21]. If p < oo, we employ 
(21), ( [22] ), ([7]) and the definition of £j to rewrite our elements: 

V0(o) • tj + \ < D 2 m^ 6 > 



= a /i I V0(O) I cos 0j + a 2 /i 2 <9i 2 0(0) cos 0, sin 7 - 

+ J a 2 h 2 A^0(O) cos 2 % + - a 2 h 2 Af 0(0) sin 2 
2 2 J 

= ah I V0(O) I cos 9j + a 2 /i 2 <9i 2 0(0) cos 9j sin % 

+ \a 2 h 2 (A«0(O) - Af 0(0)) cos 2 0, + ^a 2 /> 2 A? 0(0) 

= a/i I V0(O) I cos 9j + a 2 h 2 d 12 0(0) cos 9j sin 9j 

+ \a 2 h 2 (A« 0(0) - Af 0(0)) (cos 2 9 3 - I) + i« 2 /> 2 A p G 0(O). 



In this way, due to Lemma 7.1, we will prove our conclusion if we show 
that when h and A9 tend to zero, the p-average of 

h a M 2,2 ^12 0(0) fl . fl , 2 , 2 A^0(O) - Af 0(0) 
a/i cos^- + a /i | V ^| cos ^ sm9 3 + a h 2 | V 0(O)| C ° S j ' 

2 A£0(O) - Af 0(0) 1 

times -^-rr; , tends to — , . . , . 

a 2 h 2 ' |V0(O)| q 
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By definition of p- average, we then need to compute the argmin of the 
function Z(t): 



j=l..M 



r, 2.2^12 0(0) Q . Q 

a h cos Uj + a h cos 9j sin Uj 



+ a 2 h 2 



Ag 0(0) - Af 0(0) 
2 |V0(O)| 



|V0(O)| 
■ cos 2 9j — t 



but: argmin Z(t) = a 2 /i 2 argmin z(t), with 



it) = A9 

j=l..M 



cos 9j + a/;. 



012 0(0) fl . fl 
,„ , 7-t , cosfe 1 ,- sin0,- 
|V0(O)| J J 



+ a ft 



A^0(O) - Af 0( O ; 2 
2 |V0(O)| 



cos 2 9j — aht 



(24) 



We set a = d 12 0(O)/|V0(O)| and 6 = (A^0(O) - A^0(O))/(2 |V0(O)|) 



and recall equation (46) from the Appendix to derive: 



M 



z'(t) = —p A9 a h [(cos 9j + aha cos 9j sin 9j + ah b cos 2 9j — ah t) 

3=1 

■ | cos 9j + aha cos 9 3 - sin 9j + ahb cos 2 9j — a h t | p ~ 2 ] ; 



we next use (47) (which holds in the classical sense if p > 2 and in the weak 



sense if 1 < p < 2) and the fundamental theorem of calculus to see that for 
any d and e: 



-p \d + e\ p - 2 (d + e) = -p \d\ p ~ 2 d-pip-1) [ + 

Jd 



r 2 d S , 



hence 



z\t) = ahA9 £ [-P I cos % cos #i 

i=l..Af 

/•cos 9j + aha cos ^ sin 9j + ahb cos 2 6>j — aht 
-p{p-l) I ' \s\ p - 2 ds}. 

J COS 9n 
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And the change of variable s = cos9j + ahu gives 



z'(t) = ahA9 [-J5|cos^| p - 2 cos^- (25) 

j=l..M 

i-a cos 9j sin 9j + b cos 2 9j — t 
—ahp{jp — l) / \ cos 9j + ahu\ p ~ 2 du}. 





We remarked previously that for our original choice of m, if a direction rj 
is in the set of admissible directions so is its opposite, its orthogonal, as well 
as its reflections in each of the axes. After rotating the coordinate system, 
so that V0(x) = |V0(O)| ex, we can not make this claim any longer, as for 
example we lose the reflections with respect to the axes. Nevertheless, given 
9j, we still have 9j + tt, so that the sum of the first term in (25) is zero, and 
we obtain 

rd cos 9j sin 9j + b cos 2 9j — t 
z'(t) — —ahp(p — 1) / \cos9j + ahu\ p ~ 2 du. (26) 

J 

The argmin of z(t), call it to, is bounded by a constant independent of 
ah. This can be seen by noticing that aht is the p-average of the values 
{cos 6^- + ahacos9j sin9j + ahb cos 2 9j}; but for our choice of angles the 
p-average of the values {cos^} is zero (values are symmetric about zero), 
hence Lemma 7.3 implies |a/i£o| < Q;/i|o| +ah\b\, that is to < \a\ + |6|. 

We would like to show that the to is equal to - up to an order of 0(A9) + 

Q 

0(e). To prove our claim, we set C(a,b,q) = 2 max <\a\ + \b\, — 1, so that 

b q 
both to and - belong to the interval \t\ < C(a,b,q), and the upper limit of 
Q 

integration in (26) verifies \acos9j sin9j + b cos 2 9j — t\ < \a\ + \b\ + C(a, b, q) 
for \t\ < C(a, b, q). 

On the other hand, when p > 2, by uniform continuity if \u\ < \a\ + \b\ + 
C(a, b, q), for any e > there is a 5 t := 6 e (a, b, q) such that for < ah < 5 t 
it holds | cos 9j + ahu\ p ~ 2 = \ cos^| p " 2 + 0(e). 

Therefore, as long as \t\ < C(a,b,q), for a fixed e > 0, there is a 5 t such 
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that if < ah < 5 e , we have 

m cos 9j sin 9 a + b cos 2 9j — t 
z '(t) = -a 2 h 2 A9p(p- 1) ^ / [|cos^| p - 2 + 0(e)] rfs 

„_i »/f 



j=l..M 



9-2 



-a 2 h 2 A9 p(p-l) ( ac °s9j sin^)|cos^| p - 

j=l..M 

2 /i 2 A0pG»-l) (6cos 2 ^- t)|cos^| p - 2 + a 2 /i 2 0(e) ; 

j=l..M 



-a 2i 



here we used the fact that A0^\ =1 M = 27r. 

Given an angle 0j, as mentioned before, we can not assume we still have 
the angle 9j + | as well; nevertheless, we have among our angles an approxi- 
mation of it up to order A9. In other words, given a certain direction, in our 
pool of directions we also have its reflection up to an error of order 0(A9), 
so that 

z'(t) = a 2 h 2 0(A9) + a 2 h 2 0(e) 

- a 2 h 2 A9p{p- 1) ( b cos20 j ~ t)\ cos 6j\ p - 2 . (27) 

j=l..M 

To proceed in our proof, we recall the elementary equality 



(cos0) p d9 



i 1 1 1 

for any 1 < p < oo, and — I — = 1, 



J 2 (cos 0) p ~ 2 d0 ^ ^ ^ 



which implies 



*'(-) = a 2 /i 2 O(A0) + a 2 h 2 0(e), (28) 



since from it we deduce 



A9 ( b cos20 j - 'J I cos^| p - 2 = 0(A9). 

7=1. M ^ ^' 
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Next we notice that for any c > with \b/q ± c o (O(A0) + 0(e)) | < 
C(a,b,q), we can use (27) and (28) to obtain: 



z! (-±c (O(A0) + O(e)) 



a" 



J ±c (O(A6) + O(e))a 2 h 2 A0p(p- 1) \ cos9 j\ P 
h 2 (0(AO) + 0(e)) j l±cbA0p(p-l) ^ |cos^| p - 2 ) 

\ i=l..M / 



(29) 



We can also find a A , such that for any AO < A , there is a positive 
constant c% independent of AO and e, for which 

Cl < AO p(p-l) Icos^f" 2 . 

j=l..M 



Hence, since for any AO and e small enough to have \0(AO) + 0(e) | < 



Cl 



C(a,b,q) 



\b\ 



1 2 

we can pick a Co for which — < cq < — and \b/q ± 

Ci Ci 



c (O(A#) + 0(e)) | < C(a,b,q), from (29) we obtain 



c |O(A0) + 0(e)| < , z'(- + co \0(AO) + 0(e) \ > . 



Recalling that to = argmin z(t) is its only critical point, and z'(t) > for 



every t > to, while z'(t) < for t < t 0} see Remark 7.1 we conclude that 



argmin z(t) 



< -\O(A0) + O(e) 



(30) 



for every AO and e small enough. 

Therefore, we have that for e and AO small enough, there is a <5 e such that 
if < ah < 5 f then 



2 . Ag^O) - Af^(O) 1 
— argmin Z(t) 



< \0(A6) + 0(c) 



and the theorem follows for the case V0(x) ^ 0. 
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Assume V0(x 

r 2 

p-average of 



0, then by Definition 2.2 we need to show that the 



a 2 h 2 \ 2 



1 



< D>(0)fj, & > tends to A£ 0(0) as /i and A9 tend 



to 0. By our choice of angles, £j = ah (cos jA6, sin jA9), so that 



a 2 /i 2 V 2 



< £>V(0)&,& > = <9n0(O) cos 2 (jA0) 



+<9 22 0(O) sm 2 (jA9) + 2 <9i 2 0(O) cos(jA#) sin(jAfl) 
= l - A 2 0(O) + ^(0ii0(O) - d 22 0(O)) cos(2j A0) 
+2 d 12 0(O) cos(jA0) sin(j'A0), 

but by our assumptions if 9 is in our selection of angles, so is 9 + |, thus the 
p-average of 

^(<9n0(O) - <9 22 0(O)) cos(2jA#) + 2 9i 2 0(O) cos(jA9) sm(jA9) 



is zero, being a set of symmetric data with respect to 0. Therefore, using 
2 



@9b, 



A 



a 2 h 2 \2 



1 



< ^0(0)^,0 > 



-AoMO) 



A 2 ? 0(O) 



and the theorem is proven. 



□ 



We are now ready to show consistency of our approximation scheme. 
Before doing so we remark that although the definition of consistency we use 
is slightly different from the one given in [2], their convergent result applies 
also for this formulation. 

Theorem 4.4. Let p > 2. Our approximation scheme is consistent, that is 
for all x G Q and G C°°(Q), we have that 

limsup5(p,y,0(y)+^,0 + O <hmsupG p ( J D 2 0(y), J D0(y),0(y),y) (31) 

p— >0 y— >x 
y ^ x yen 

(where, as before, p = (h,A9)) and 
liminf5(p,y,0(y)+^,0 + O >hminfG p ( J D 2 0(y), J D0(y),0(y),y). (32) 



p->0 

y— >x 



y-»c 
yen 
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Proof. If x G Q both statements are a consequence of our previous theorem, 
i.e. Theorem 4^3 On the other hand, if x G dfl, we have that 

hmsupG p ( J D 2 0(y), J D0(y),0(y),y) 



yen 



max (H p (D 2 0(x) , Zty(x) , x) , 0(x) - F(x)) 



while 



liminfG p (DV(y), J D0(y),0(y),y) 

yen 



= min (^(^(x),^(x),x)^(x) - F(x)) , 

and the theorem follows by the definition of 5". □ 

However, as far as an explicit scheme in time is applied to 

u t + G p (D 2 u(y), Du(y), u(y), y) = 0, (33) 

the consistency of the scheme S with respect to the stationary nonlinear 
operator G p implies the consistency for the evolutive operator as in [2]. 
In fact, by applying the Euler approximation in time we have, for a given 
initial condition u°, the explicit time marching scheme 

u n+l =u n_ At S fa y> U »(y), U ») } (34) 

which implies, taking At = \p\, 



u n+l _ u r 



\P\ 



-S(p,y,u n (y),u n ). 



Passing to the limit for \p\ which tends to 0, we get consistency in the usual 
sense. 

Theorem 4.5. Let f = 0. For all h > 0,A6 > 0, there exists a solution 
Up e L°°(f2) of (15) such that \ \u p \\ L oo(n-\ < | |z,<»(an) ■ 

Proof. We consider the operator E p : L°°(f2)) — > L°°(0)) defined as 

A p {C£ e {x,u;a)) if x G Q 

F(x) if x G dn, 
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and notice that thanks to Lemma 



7.1 



in the Appendix, we have 



\E p u\ 



L°°(Q) 



< max d I it I 



L°°(n)> 



|L°°(dfi), 



(35) 



Additionally, E p is a nonexpansive operator in the L°° norm, that is 



\E p (u)-E p (v) 



Il°°(q) 



li°°(n)) 



(36) 



since if x e f2, by Lemma 7.3 we know that 

\E p (u)(x)-E p (v)(x)\ = \A p (Ct e (^u;a))-A p (Ct e (^v;a))\ < \\u-v\ 



L°°(n)' 



while 



\E p (u)(x) - E p (v)(x)\ =0, if x e dQ. 



If we consider B(0,\\F\\ L oo^ 9n ^) to be the sphere centered at the zero 
function and of radius H-FH^^n), equations (35) and (36) imply that E p is a 
nonexpansive operator mapping the closed sphere B(0, \ \F\\L°°(dn)) C L°°(Q) 
into itself. Therefore, by a classical fixed point theorem result (see Corollary 1 
and Remark thereafter in [26]) we conclude that E p has a fixed point in 
B(0, \ \F\\ L oo(Qfy), and the theorem follows. □ 

Theorem 4.6. Assume f = and p > 2, then the solution u p of the appro- 
ximation scheme (15) converges as p = (h, AO) — > (0,0) to the viscosity 
solution u of (19). 



Proof. For / = 0, we know that (19) has a unique bounded viscosity solution 
u, with IMIioo/Q) < ||^||i«>(9n) and that a comparison principle holds (see 



[T9l [16] for details). Therefore, thanks to Theorems 4.2, 4.4 and 4.5 
apply Theorem 2.1 in (2]. 



we can 

□ 



5. Numerical implementation 

Let us now introduce in Q a structured grid according to a space dis- 
cretization parameter h, with < h < 1, and denote by {x.j}j=i..N its nodes. 
An important step to implement our approximation scheme for a fixed angu- 
lar resolution AO is the way we reconstruct the values (x, it; a) in (14), 
starting from the known values of u at the nodes of the grid. This is done 
via interpolation, but in order to obtain a convergence result we must re- 
strict to monotone interpolation techniques. This means that denoted by u 
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a function defined in a domain D and by /[«](•) its local interpolation based 
on the values at the nodes, we can only use an interpolation operator such 
that 

m = minw(:r) < < maxu(x) = M. (37) 

Note that linear and bilinear interpolation in IR 2 are monotone interpolation. 
Another important property of bilinear interpolation is that it is translation 
invariant, i.e. given a constant S we have 

I[u + 5](x) = I[u] + 5. (38) 

This property will guarantee that the resulting scheme is also invariant with 
respect to the addition of constants. 

We implement our approximation scheme on a structured grid using the 



classical time marching approximation (34). Since our solution u n is charac- 
terized by its values at the nodes, we prefer to work in the space M. N , denoting 
by u n the iterate at step n, that is the vector of components = u n (xj); 
given an initial condition u° e R N we then consider the iterative scheme 



2 At 

n V> 2 



A p (d^(x,-,u n ;%))- < +A*/( Xi ) x 3 GH 



u] +1 =\ i (39) 

F(xj) Xj G dtt; 

if for a given integer m we assume A6 = ir/(2m), now 

a^(x„ u"; a 3 ) = { W (x}, u n ),i = 0, 4m - 1}, (40) 

where x* = Xj + hajTi (with r-j = (cos iA8, sin i AO)), and w(xj,u n ) is the 
value of the function u n at x* obtained by using bilinear interpolation on 
the four closest grid-points (in fact this is the main difference with respect 
to C^ e ). Here, aj is a parameter that may vary at each grid point x^, and 
which verifies < a* < ctj < distfej, dQ) < dn, for some constant a* and 
for dn, the diameter of Q. In other words, acj is a function of j, which is 
uniformly bounded independently of h from above and below, and such that 
x*- belongs to the computational domain whenever < h < 1. 

More precisely, if x = (x, y) is a point contained in a grid cell with vertices 
(x k ,yi) (the lower left corner), (x k+1 ,yi), (x k ,y l+1 ) and (x k+1 ,y i+1 ), and we 
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denote by u k) i = u(x k ,yi), etc., the known values of a function u on them, 
the bilinear interpolation computes in (x, y) the second order polynomial 

I[u](x, y) = axy + bx + cy + d, 

where the coefficients a, b, c and d can be determined solving the linear 
system 4x4 which corresponds to the four height conditions at the four 
vertices of the cell. Those values can also be written as linear combinations 
of the values of u at the vertices of the cell, i.e. 

I[u)(x,y) = \k,iUk,i + ^k+i,iUk+i,i + ^k,i+\Uk,i+i + ^k+i,i+i u k+i,i+i, 

where the coefficients are given by 

A fc ,z = (x fc+ i - x)(y i+1 - y)/V, 

Afc+i,/ = (x- x k )(y l+1 - y)/V, 

= (xfc+i -x)(y- yi)/V, 

\ k +i,i+i = (x- x k )(y - yi)/V, 

and V := (x k +i — X)-)(yi+i — yi) is the area of the cell, i.e. V = h 2 for our 
uniform grid (see for example [T3]). 

The p-average in (39) is then computed through the 4m values of I[u n ] 
on a set of equally distributed points on the sphere of center x^ and radius 
hazj. 

Then the approximate solution of the Dirichlet problem ^ for the game 
p-Laplacian is computed by using the scheme above and running it until the 
stopping rule 

E n = max|M™ +1 -u n A < e (41) 



is satisfied, for a given tolerance e. Another option is to use the simpler 

2 At r~i 

iterative scheme obtained by setting 2 = 1 in (39), that is 



a]h 2 



ur 1 - 



if Xj G Cl, 



if Xj G dil; 



(42) 



until convergence (i.e. until (41) is satisfied). 
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We next show that for a suitable initial configuration the iteration gener- 
ated by (39) converges if / = 0. Although we can not claim that this proves 
convergence to the solution of our approximation scheme, the numerical tests 
that we present in Section [6] show convergence to the correct viscosity solu- 
tion in all cases were the exact solution is known, also for / ^ 0. Let us 
rewrite our numerical iteration as 



u 



.71+1 



where for u G 



we have set 



Uj + 



2 At 



(43) 



u , 



+ Atf{x j ), x 3 GH, 



(44) 



Xj e dCl. 



It is not difficult to show that for a given grid, the fact that the boundary 
nodes have a fixed value over each iteration prevents the numerical solution 
from blowing up even in the presence of the source term /. On the other 
hand, for / ^ the bound depends on the grid size and / jointly, in a 
manner for which we don't have an independent bound. Instead, if / = it 
is very easy to derive that the initial condition u° will provide a bound for 
any subsequent iteration as the following result shows. 
2 At 

Theorem 5.1. Let -^r^ < 1. Assume f = 0, then for n > 1 it holds 



a 2 h 2 
sup 

j=l..N 



< sup 

j=l..N 



U 



n-l\ 



< sup 

j=l..N 



II : 



where u" is defined by (43). 

Proof. It is enough to look at the internal nodes, since the approximation 
at the boundary nodes has fixed values. For j fixed, the set Cfc (x 
consists of values w(x 
Hence, each w(jcI 



3 1 



u 11 ) computed by using bilinear interpolation of u n . 

and Lemma 



-j,u") is controlled by the values u^, 
Appendix implies that 



7.1 



in the 



< 1- 



\( T p( u )) 



2 At 
1 ~ ^h? 
2 At 



u 



n-l 



\U 



n-l i 



2 At 
afh 2 
2 At 



sup \u 

j..N 



n-l I 
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and stability follows. □ 

The next theorem shows that if the initial condition is appropriately cho- 
sen, then the iterates generated by our schemes are point-wise increasing; this 
fact together with the previous stability result implies that they pointwise 
converge if / = 0. 

2 At 

Theorem 5.2. Let < 1. For n > 1, there exists an initial condition, 

a% h 2 

u°, for which the iterations generated by the scheme (43) verify 

> f or an V j = 1— ^V> an d n > 1. 
Proof. We choose as initial condition: 



mini* 1 if x,- G f2, 

an J 



and since vR > min F for any j we have 



an 

MCffo, u°; > min F if x, G fl. 

oil 



Therefore, 



u) ={l- u] + ^Mdfte, u°; aj )) + At f( Xj ) 

> min F + At fixA > u° H , if x,- G Q, 
since / > 0. Then, we conclude it] > w° for every j, as well as 

^(^(x,, u 1 ; a 3 )) > MCfixj, u°; a,)), if x, G fl, 



thanks to Lemma 7.2 Inserting this last inequality in the definition of u 
leads to 



2 At 

a 



u 2 = { 1 - «} + ^^(^(x^u 1 ;^)) + At f(x 



w], if Xj G f2, 



23 



that is v% > vh for every j. We then obtain the desired conclusion by 
induction. □ 

Remark 5.1. In the case f = 0, if we pick as initial condition 



4 



maxF if x,- G £X 

an 



and follow the steps of the proof of Theorem we obtain a non-increasing 



sequence, which dominates element by element the sequence in Theorem 5.2 



Note that, if these two sequences converge to the same limit u, choosing 
as initial iteration a u° such that u® is between the minimum and maximum 
values of F, we would have a sequence converging again to u. Since for f = 
we know that problem (jlp has a unique viscosity solution, if we could show 
that our numerical implementation converges to it then any initial condition 
such that mingn F < < maxgn F would produce a sequence converging 
to the viscosity solution. Let us also observe that the choice of a monotone 
interpolation in the numerical implementation guarantees that similar bounds 
also apply to the interpolation of our initial condition and to all the elements 
in the sequence. 

6. Numerical tests 

We present in this section some experiments obtained with our numer- 
ical implementation coded in MATLAB, and executed on a MacBook Pro 
desktop machine with a 2.2 GHz Intel Core 2 Duo processor. As described 
in Section 5, we have used structured uniform grids, and the values of the 



approximate solution at the points in C^ e of (40) have been computed via 
bilinear interpolation using the four closest grid points. If one of the points 
lies on a line joining two grid points, the bilinear interpolation reduces to 
the linear interpolation between them. To compute the p-average at each 
node we used the Newton Bracketing method for minimization of convex 
functions, [IS] , applied to the function g(s) := (Q(s, S))?, where Q(s,S) is 



defined as in equation (45). We believe that an optimization of this part of 
the procedure could improve the speed of calculations. 

Depending on the examples we have picked different values for the param- 
eter a j at different grid points. To better illustrate our choices, we introduce 
the following simple definition, where dj denotes the distance of xj from dQ 
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Figure 1: The 4-level circles stencil, with j3 = 1 



along the grid lines, and therefore -j- is an integer greater or equal to 1 (see 
Figure 1). 



Definition 6.1. An iteration generated by (39) is called n-level circles iter- 
ation if for every j the parameter aj is chosen so that aj = (3mm(n, -jr), for 
a given < (3 < 1. 



Test 1. For comparison, first of all we run the schemes of the previous 
section on the examples for the oo-Laplacian presented in [2T|, using the 
same specifications provided there (in this case convergence is assumed to be 
reached when E n < 2/il0~ 2 , where E n is the quantity defined in (41)). 
We have applied scheme (39) to the case Q = (—1, 1) x (- 



-1,1), P 



oo, 

/ = 0,F(x,y) = \x\ 4 / 3 — and the results are summarized in Table 1 

below. The explicit solution of this problem is the well-known Aronsson 
function u(x,y) = \x\ 4 ^ 3 — \y\ 4 ^ 3 . We denote by iV 2 (so that h = 2/N) the 
total number of nodes in the square grid, and by n the number of iterations 
(reported between parentheses), while the error is given in the maximum 
norm and a 2-level circles iteration is used, for = 0.99. The initial con- 
dition u° was assumed to be a perturbation of the exact solution (of order 
±20%). The total number of grid points and the number of directions used 
to compute Cfr 9 have been choosen to make our tests comparable with the 



ones for the 17-point and the 25-point stencils presented in Table 2 of [2U 
pg. 1227]. Note that the errors decrease on the rows and on the columns in 
a regular way and that also the number of iterations decreases if we increase 
the number of directions and iV simultaneously. Our numerical results show 
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essentially the same accuracy of those of [21J. We only remark that in our 
case, in order to significantly reduce the error when h tend to zero, we should 
also increase the number of directions. 



Dir. 


iV = 41 


N = 81 


N = 161 


N = 241 


N = 401 


4 


0.1105 (250) 


0.0765 (448) 


0.0373 (584) 


0.0225 (589) 


0.0122 (621) 


8 


0.0274 (80) 


0.0182 (161) 


0.0084 (214) 


0.0069 (190) 


0.0048 (188) 


16 


0.0084 (54) 


0.0070 (75) 


0.0043 (105) 


0.0033 (108) 


0.0023 (112) 


24 


0.0088 (57) 


0.0081 (73) 


0.0050 (91) 


0.0035 (103) 


0.0024 (107) 



Table 1: L°°- errors and iterations (in parentheses) for Test 1 on Aronsson function 



In Figure 2 one can see the numerical solutions and its contour plots 
obtained, when again Q = (—1,1) x (—1,1), p = oo and / = 0, for three 
different choices of the boundary data. These computations are included for 
comparison with Figure 2 in [2U pg. 1228]. In the first two cases, we consider 
F(x,y) = \x\ 2 \y\ 2 and F(x,y) = x 3 — 3xy 2 . We use 24 controls, 2- level circles 
iterations, /3 = 0.99, and a grid with 401 2 nodes. The initial condition for 
this grid is generated by a multi-resolution type approach. More specifically, 
we start from a 21 by 21 coarse grid with initial values set to zero in the 
interior nodes; after a few iterations we interpolate the numerical solution on 
a finer grid and repeat the procedure up to the desired resolution. Although 
this start up procedure requires some time, but significantly improves the 
rate of convergence. The approximations in Figure 2 for these two cases 
took n — 15 and n = 47 iterations, respectively, to converge. In the third 
example, F is the characteristic function of the point (1,0). We proceed as 
described above, but we use 24 controls on a J^-level circles iteration, again 
for (3 = 0.99. In this test, the approximation took n = 2428 iterations to 
converge with the required accuracy. 

Test 2. We next consider the problem for the game p-Laplacian in a case 
where we were able to compute the exact solution. Starting with Q = B(0, 1), 
/ = 1 and F = 0, and working in radial coordinates we derived the solution 

i 2 2 

for any p > 2, that is v(x,y) = ~ x 2 ~ y ■ This is an analytic function in 
the whole plane and has a unique extrema at the origin. By looking at v 
on the unit square Q = (—1,1) x (—1, 1), we see that this function verifies, 
for any p > 2, — A^v = 1 in Q. To test our code, we have implemented it 

on the problem: — A^u — 1 in Q, for F(x, y) = 1 ~ x 2 ~ v on the boundary. 
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Figure 2: Test 1, numerical solutions and corresponding contour plots when Q = ( — 1, l) 2 , 
p = oo, / = 0, with boundary data F(x,y) — \x\ 2 |y| 2 , F(x,y) — x 3 — 3xy 2 , and F the 
characteristic function of point (1,0). 
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We summarize the numerical results obtained with scheme (42) for the cases 
p = 5 and p = oo in Tables 2 and 3 below, showing the L°°-errors and the 
number of iterations of the algorithm until convergence (E n < 10~ 5 , for this 
test) for different combinations of levels and directions. The initial iteration 
was set to minF = —1/2 at the interior nodes. 



Nodes (levels) 


16 directions error (iter) 


24 directions error (iter) 


21 (2) 


0.0634 (163) 


0.0617 (180) 


21 (4) 


0.0241 (50) 


0.0192 (107) 


41 (4) 


0.0201 (213) 


0.0191 (163) 



Table 2: Test 2, L c 
(3 = 0.9. 



-errors and iterations for p — 5, / = 1, F(x,y) = (1 



2/ 2 )/2, 



Nodes (levels) 


16 directions error (iter) 


24 directions error (iter) 


21 (2) 


0.0590 (249) 


0.0563 (248) 


21 (4) 


0.0211 (80) 


0.0185 (77) 


41 (4) 


0.0192 (272) 


0.0156 (272) 



Table 3: Test 2, L°°-errors and iterations for p = oo, / = 1, F(x,y) = (1 — x 2 — y 2 )/2, 
(3 = 0.9. 



Test 3. ( The tug-of-war game) We finally run the code on the rectangle 
= (—2,2) x (—1, 1), for different values of p, with / = 1 and F = 0. In 
particular, the case p = oo and / = 1 corresponds to a running cost in the 
tug-of-war game. The exact solution is not known, and to our knowledge 
these are the first numerical results for a solution without radial symmetry. 
More precisely, the explicit solution is known only in a part of the domain, 
where it can be computed based on ideas from [21]. That is, for any — 1 < 
x < 1, the exact solution is given by u{x, y) = We recover such solution 
in this part of the domain with our numerics. We set the initial condition 
u° to minF = in the interior nodes, and run the scheme M2J) with J^-level 
circles iterations, with (3 = 0.8, until the difference of the values at (0,0) of 
two consecutive approximations is less than 10~ 6 (note that u(0, 0) = 0.5 is 
the maximum value of the exact solution). Our tests have shown that the 
scheme converges for any choice of interior values for u°, with of course only a 
different number of iterations. Here we summarize in Table 4 and in Figure 3 
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below the results of the simulation with 16 controls (due to the axis-oriented 
solution, there is no real advantage in this case to use more directions). 

We have also included Table 5 in order to compare the multi-level circle 
approaches: as expected the table shows some acceleration of convergence 
with larger multi-level circles, since the information from the boundary can 
reach the interior of the domain quicker for larger circles. 




grid 


error at (0,0) 


iterations 


CPU time 


161 x 81 


0.0276 


1112 


236 


241 x 121 


0.0155 


2205 


957 


321 x 161 


0.0094 


3578 


2681 



Table 4: Test 3, p = oo, fi = (-2, 2) x (-1, 1), / = 1, F = 0, /3 = 0.8, 16 directions. 



7. Appendix: some p-average properties 

For reader's convenience we provide here some elementary properties sat- 
isfied by the p-average of finite sets of real numbers. 

For a fixed set S = {si,...,s m } C E, s e K and p > 1, we define the 
function 

m 

Q(s,S) = J2\ s i- s \ P > ( 45 ) 

3=1 
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levels 


error at (0,0) 


iterations 


CPU time 


4 


0.0276 


1112 


236 


2 


0.0260 


3330 


621 


1 


0.0917 


9206 


1881 



Table 5: Test 3, p = oo, ft = (-2,2) x (-1,1), / = 1, F = 0, $ = 0.8, 16 directions, 
161 x 81 nodes. 



whose derivative with respect to s is easily computed as 
dQ, 



ds 



(a, S) = - P J2 \'i ~ sgn( S , - s) = p £ \s - stf' 2 (s - 8j ). (46) 



For any p > 1 one can also compute the second derivative of Q(s, S) with 
respect to s: 

^(s,S)=p(p-l)J2\s 3 -sr 2 . (47) 

In the case 1 < p < 2 such relation has to be understood in the weak sense 
of W 2 ' 1 functions. 

Remark 7.1. The function Q has exactly one extremum and is convex in 

dQ 

s, and the p-average A p (S) is the only value for which -^-(A p (S), S) = 0, 

dQ dQ 
hence -^-(s,S) < for every s < A P (S), and -^-(s,S) > for s > A P (S). 

Additionally, A P (S) solves implicitly the following equation: 



p-2 



A (S) = ^ 8 ^Ms)\°3 ^pw;r Sj 

Y. S:j ^A p (s) \ s j ~ A p{s)\ p ~ 2 

Lemma 7.1. Let S = {si, S2, s m } be a finite set of real numbers, and for 
I: el let S + k = {si + k, S2 + k, s m + k}. The following assertions hold 
true for 1 < p < oo: 

A P (S + k) = A P (S) + k, (49) 
min Sj < A p (S) < max Sj . (50) 

j=l..m j=l..m 
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Proof. The cases p = 1, p = oo are trivial. For p > 1, the first assertion 



follows by the uniqueness of A P (S), while the second follows by Remark 7.1 
and the fact that if s„ 



min 

j=l..m 



and s* 



max Sj then by (|46|) one has 

i=l..m 



^-(s*,S) <0and ^(s*,5) > 0. 
as as 



□ 



Lemma 7.2. Let S = {si, s 2 , s m } and T = {ti,t 2 , ...,t m } be two finite 
sets of real numbers having the same number m of elements, and let 1 < p < 
oo be fixed. If it holds that tj 
A P {T)<A P {S). 



< Sj, for every j 



m, then we have 



Proof. According to the definition of p-average given in (11) the lemma 
clearly holds for the cases p — 1, and p = oo. 



Assume 1 < p < oo. 



For r G 1Z define the function M p (r) 



'r if 



r / and M p (0) = 0, note that MJr) is a continuous increasing function. 



Given t := A p (T), by Remark 7.1 we know that ^P-(t,T) = 0, thus by 



equation (46) we obtain Y2T=i M p (i — tj) = 0. But 



ds 

since t~ 



< Sj it holds 



< 



M p (t - tj), therefore Z?=i M p (t - Sj ) < £7=i M p (t - t 



M p (t- 

which gives f£(t, S) < 0. By Remark fu\ we conclude A p (S) > t := A p (T). 

□ 



Lemma 7.3. Le^ S 1 and T be two finite sets of real numbers having the 
same number of elements, and let 1 < p < oo be fixed. Assume that S = 
{s 1 ,s 2 ,...,s m } and T = {t h t 2 , t m } verify tj = Sj + 8j, for every j = 
1, ...,m, where \5j\ < 5 for some 5 > 0, then one has 



Ap(S) -S< A P {T) < A P {S) + 5. 



(51) 



5 < tj < Sj + 5, Lemma 7.2 implies A p (S — 5) < A P (T) < 



Proof. Since Sj 
A p (S + 5). But, from equation (49) we know A P (S — 5) 
A p (S + 5) = A P (S) + 5, so that fan follows. 



A p (S) — 5, and 
□ 
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